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The Metric Traveling Salesman Problem (TSP) is a classical NP-hard optimization problem. The 
double-tree shortcutting method for Metric TSP yields an exponentially-sized space of TSP tours, 
each of which approximates the optimal solution within at most a factor of 2. We consider 
the problem of finding among these tours the one that gives the closest approximation, i.e. the 
minimum-weight double-tree shortcutting. Burkard et al. gave an algorithm for this problem, run- 
ning in time 0(n^ + 2'^n^) and memory 0(2'^n^), where d is the maximum node degree in the 
rooted minimum spanning tree. We give an improved algorithm for the case of small d (including 
planar Euclidean TSP, where d < A), running in time 0(4'*n^) and memory 0(4'^n). This improve- 
ment allows one to solve the problem on much larger instances than previously attempted. Our 
computational experiments suggest that in terms of the time-quality tradeoff, the minimum- weight 
double-tree shortcutting method provides one of the best known tour-constructing heuristics. 

Categories and Subject Descriptors: E.l [DATA STRUCTURES]: Graphs and networks; F.2.2 
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MATICS]; Graph Theory — Graph algorithms 
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1. INTRODUCTION 

The IVIetric Travelling Salesman Problem (TSP) is a classical combinatorial opti- 
mization problem. We represent a set of n points in a metric space by a complete 
weighted graph on n nodes, where the weight of an edge is defined by the distance 
between the corresponding points. The objective of JVIetric TSP is to find in this 
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graph a minimum-weight Hamiltonian cycle (equivalently, a minimum-weight tour 
visiting every node at least once). The most common example of Metric TSP is 
the planar Euclidean TSP, where the points lie in the two-dimensional Euclidean 
plane, and the distances are measured according to the Euclidean metric. 

Metric TSP, even restricted to planar Euclidean TSP, is well-known to be NP- 
hard [Papadimitriou 1977]. Metric TSP is also known to be NP-hard to approximate 
to within a ratio 1.00456, but polynomial-time approximable to within a ratio 1.5. 
Fixed-dimension Euclidean TSP is known to have a PTAS (i.e. a family of algo- 
rithms with approximation ratio arbitrarily close to 1) [Arora 1998]; this generalises 
to any metric defined by a fixed-dimension Minkowski vector norm. 

Two simple approaches, the double-tree method [Rosenkrantz et al. 1977] and the 
Christofides method [Christofides 1976; Serdyukov 1978], allow one to approximate 
the solution of Metric TSP within a factor of 2 and 1.5, respectively. Both methods 
belong to the class of tour-constructing heuristics, i.e. "heuristics that incremen- 
tally construct a tour and stop as soon as a valid tour is created" [Johnson and 
McGeoch 2002]. In both methods, we build an Eulerian graph on the given point 
set, select an Euler tour of the graph, and then perform shortcutting on this tour 
by removing repeated nodes, until all node repetitions are removed. In general, it 
is not prescribed which one of several occurrences of a particular node to remove. 
Therefore, the methods yield an exponentially-sized space of TSP tours (shortcut- 
tings of a specific Euler tour in a specific Eulerian graph), each approximating the 
optimal solution within a factor of 2 (respectively, 1.5). 

The two methods differ in the way the initial weighted Eulerian graph is con- 
structed. Both start by finding the graph's minimum- weight spanning tree (MST). 
The double-tree method then doubles every edge in the MST, while the Christofides 
method adds to the MST a minimum- weight matching built on the set of odd-degree 
nodes. The weight of the resulting Euler tour exceeds the weight of the optimal TSP 
tour by at most a factor of 2 (respectively, 1.5), and the subsequent shortcutting 
can only decrease the tour weight. 

While any tour obtained by shortcutting of the original Euler tour approximates 
the optimal solution within the specified factor, clearly, it is still desirable to find the 
shortcutting that gives the closest approximation. Given an Eulerian graph on a set 
of points, we will consider its minimum-weight shortcutting across all shortcuttings 
of all possible Euler tours of the graph. We shall correspondingly speak about the 
minimum-weight double-tree and the minimum-weight Christofides methods. 

Unfortunately, for general Metric TSP, both the double-tree and Christofides 
minimum-weight shortcutting problems arc NP-hard. Consider an instance of the 
Hamiltonian cycle problem on an unweighted graph; this can be regarded as an 
instance of Metric TSP with weights 1 and 2. Add an extra node connected to 
all the original nodc;s by edges of weight 1, and take the newly added edges as the 
MST. It is easy to see that the resulting minimum-weight double-tree shortcutting 
problem is equivalent to the original Hamiltonian cycle problem. The minimum- 
weight double-tree shortcutting problem was believed for a long time to be NP-hard 
even for planar Euclidean TSP, until a polynomial-time algorithm was given by 
Burkard et al. [1998]. This is the algorithm we improve upon in the current paper. 
In contrast, the minimum- weight Christofides shortcutting problem remains NP- 
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hard even for planar Euclidean TSP [Papadimitriou and Vazirani 1984]. 

In the rest of this paper, we will mainly deal with the rooted MST, which is 
obtained from the MST by selecting an arbitrary node as the root. In the rooted 
MST, the terms parent, child, ancestor, descendant, sibling, leaf all have their 
standard meaning. Let d denote the maximum number of children per node in 
the rooted MST. Note that in the Euclidean plane, the maximum degree of an 
unrooted MST is at most 6. Moreover, a node can have degree equal to 6, only if it 
is surrounded by six equidistant nodes forming a regular hexagon; we can exclude 
this degenerate case from consideration by a slight perturbation of the input points. 
This leaves us with an unrooted MST of maximum degree 5. By choosing a node 
of degree less than 5 as the root, we obtain a rooted MST with d < 4. 

The minimum-weight double-tree shortcutting algorithm of [Burkard et al. 1998] 
applies to the general Metric TSP, and runs in time 0{n^ + and memory 

0{2'^n'^). In this paper, we give an improved algorithm^ for the case of small d, 
running in time 0{4'^n'^) and memory 0{4'^n). In the planar Euclidean case, both 
above algorithms run in polynomial time and memory. 

We then describe our implementation of the new algorithm, which incorporates 
a couple of additional hcTiristic improvements designed to speed up the algorithm 
and to increase its approximation quality. Computational experiments show that 
the approximation quality and running time of our implementation are among the 
best known tour-constructing heuristics. 

A preliminary version of this paper appeared as [Deineko and Tiskin 2007] . 

2. THE ALGORITHM 

2.1 Preliminaries 

Let G be a weighted graph representing the Metric TSP problem on n points. The 
double-tree method consists of the following stages: 

— construct the minimum spanning tree of G; 

— duplicate every edge of the tree, obtaining an n-node Eulerian graph; 

— select an Euler tour of the double-tree graph; 

— reduce the Euler tour to a Hamiltonian cycle by repeated shortcutting, i.e. re- 
placing a node sequence a, b, c by a, c, as long as node b appears elsewhere in the 
current tour. 

We say that a Hamiltonian cycle conforms to the doubled spanning tree, if it can 
be obtained from that tree by shortcutting one of its Euler tours. We also extend 
this definition to paths, saying that a path conforms to the tree, if it is a subpath 
of a conforming Hamiltonian cycle. 

In our minimum-weight double-tree shortcutting algorithm, we refine the bottom- 
up dynamic programming approach of [Burkard et al. 1998]. Initially, we select an 
arbitrary node r as the root of the tree. For a node u, we denote by C'{u) the set 
of all children of u, and by T{u) the node set of the maximal subtree rooted at u. 



-"^Note that Burkard ct al. [Burkard ct al. 1998] also give an 0{2''n'^) algorithm for a more gen- 
eral TSP-type problem, where the set of admissible tours is restricted by a given PQ-tree. Our 
algorithm does not improve on the algorithm of [Burkard et al. 1998] for this more general problem. 
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i.e. the set of all descendants of u (including u itself). For a set of siblings U, we 
denote by T{U) the (disjoint) union of all subtrees T{u), u G U. When U is empty, 
T{U) is also empty. 

The characteristic property of a conforming Hamiltonian cycle is as follows: for 
every node u, the cycle must contain all nodes of T{u) consecutively in some order. 
For an arbitrary node set S, wc will say that a path through the graph sweeps 
S, if it visits all nodes of S consecutively in some order. In this terminology, a 
conforming Hamiltonian cycle must, for every node u, contain a subpath sweeping 
the subtree T(u). 

In the rest of this section, we denote the metric distance between u and v by 
d{u, v). We use the symbol ttl to denote disjoint set union. For brevity, given a set 
A and an element a, we write A\£a instead of {a}, and A\a instead of A\{a}. 



2.2 Upsweep: Computing solution weight 

The algorithm proceeds by computing minimum- weight sweeping paths in progres- 
sively increasing subtrees, beginning with the leaves and finishing with the whole 
tree T{r). A similar approach is adopted in [Burkard et al. 1998], where in each 
subtree, all-pairs minimum-weight sweeping paths are computed. In contrast, our 
algorithm only computes single-source minimum- weight sweeping paths originating 
at the subtree's root. This leads to substantial savings in time and memory. 

A non-root node v & C{u) is active, if its subtree T{v) has already been processed, 
but its parent's subtree T{u) has not yet been processed. In every stage of the 
algorithm, we choose the current node u, so that all children of u (if any) are 
active. We call T{u) the current subtree. Let V C C{u), a & T{V). We denote by 
Dy(a) the weight of the shortest conforming path starting from u, sweeping the 
subtree u^BT{V), and finishing at a. 

Consider the current subtree T{u). Processing this subtree will yield the values 
-Dy(a) for all V C C(u), a. G T{V). In order to process the subtree, wc need the 
corresponding values for all subtrees rooted at the children of u. More precisely, we 
need the values D^^{a) for every child v e C{u), every subset W C C{v), and every 
destination node a G T{W). Wc do not need any explicit information on subtrees 
rooted at grandchildren and lower descendants of u. 

Given the current subtree T{u), the values Dy{a) are computed inductively for 
all sets V of children of u. The induction is on the size of the set V. The base of 
the induction is trivial: no values Dy{a) exist when F = 0. 

In the inductive step, given a set ^ C C{u), we compute the values L'yyj^(a) for 
all V G C(u) \V, a <E T{v), as follows. By the inductive hypothesis, we have the 
values Dy{a) for all a G T{V). The main part of the inductive step consists in 
computing a set of auxiliary values Dyyy{v), for all subsets W C C{v). Every such 
value represents the weight of the shortest conforming path starting from node u, 
sweeping the subtree u tt) T{V), then sweeping the subtree T(W) l±l v, and finishing 
at node v. Suppose the path exits the subtree u W T{V) at node x and enters the 
subtree T{W) l±l at node y. We have 
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T{V) 
(a) Case W = 9 



T{V) 
(b) Case W ^ 



Fig. 2: Computation of D^|^^(a), a e T{v) 



'd{u,v) iiV = ^,W = % 

Any^T(w) [d{u, y) + D-^{y)\ \iV = $,W^<h 

min,er(y) [Dl{x) + d{x, v)] if F ^ 0, = 

. min,eT(y).^eT(vy) [^'^^ {x) + d{x, y) + D^^ {y)] if F ^ 0, 1^ ^ 



(1) 

(see Figure 1). The required values Dy^{y) have been obtained previously, while 
processing subtrees T{v) for the active nodes v G C{u). Note that the computed 



auxiliary values include -Dyy^(t') = D 



V,C{v) 



(v). 



Now we can compute the values Dy^^^^{a) for all a e T{v) \ v = T{C{v)). A path 
corresponding to Dy^^(a) must sweep u ttJ T{V), and then T{v), finishing at o. 
While in T{v), the path will first sweep a (possibly single-node) subtree vkBT{W), 
finishing at v. Then, starting at v, the path will sweep the subtree v^T{W), where 
W — C{v) \ W, finishing at a. Considering every possible disjoint bipartitioning 
WtiiW = C{v), such that a G T{W), we have 



D^M= _ min _[DIM+D^ia) 

WiSW=C{v): aeT{W) 



(2) 



(see Figure 2). 

Wo now have the values Dy^^{a) for all a G T{v). The computation (l)-(2) is 
repeated for every node v G C{u) \ V. The inductive step is now completed. 

The processing of subtree T{u) terminates when all possible choices of subset V 
and node v have been exhausted. 

Eventually, the root r of the tree becomes the current node, and we process the 
complete tree T(r). This establishes the values Dg{a) for all S C C{r), a G T{S), 
which includes the values -D^^^j (a) for all a^r. The weight of the minimum- weight 
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T(a) 



Fig. 3: Computation of P^{a), a € T{V), k = 3 



conforming Hamiltonian cycle can now be determined as 




(3) 



Theorem 2.1. The upsweep algorithm computes the weight of the minimum- 
weight tree shortcutting in time (9(4''n^) and space 0{2'^n). 

Proof. In computation (1), the total number of quadruples u,v,x,y is at most 
(since for every pair x, y, the node u is determined uniquely as the lowest 
common ancestor of x, y, and the node v is determined uniquely as a child of u 
and an ancestor of y). In computation (2), the total number of triples u, v, a is also 
at most (since for every pair u, a, the node v is determined uniquely as a child 
of u and an ancestor of y). For every such quadruple or triple, the computation is 
performed at most 4'' times, corresponding to 2"^ possible choices of each of V, W. 
The cost of computation (3) is negligible. Therefore, the total time complexity of 
the algorithm is 0(4'*n^). 

Since our goal at this stage is just to compute the solution weight, at any given 
moment we only need to store the values Dy{a), where u is either an active node, or 
the current node (i.e. the node for which these values are currently being computed). 
When u corresponds to an active node, the number of possible pairs u, a is at most 
n (since node u is determined uniquely as the root of the active subtree containing 
a). When u corresponds to the current node, the niunber of possible pairs u,a is 
also at most n (since node u is fixed). For every such pair, we need to keep at most 
2** values, corresponding to 2"^ possible choices of V. The remaining space costs are 
negligible. Therefore, the total space complexity of the algorithm is 0(2'*n). □ 

2.3 Downsweep: Reconstructing full solution 

In order to reconstruct the minimum- weight Hamiltonian cycle itself, we miist keep 
all the auxiliary values Dy^^{v) obtained in the course of the upsweep computation 
for every parent-child pair u, v. We solve recursively the following problem: given 
a node u, a set V C C{u), and a node a E T(V), find the minimum- weight path 
Pyia) starting from u, sweeping subtree u W T{V), and finishing at a. To compute 
the global minimum-weight Hamiltonian cycle, it is sufficient to determine the path 
^C{r)i"-)' where r is the root of the tree, and a is the node for which the minimum 
in (3) is attained. 

For any u, V C C{u), a G T{V), consider the (not necessarily conforming or 
minimum- weight) path u = vq ^ Vi ^ V2 ^ ■ ■ ■ ^ Vk = a, joining nodes u and 
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a in the tree (see Figure 3). The conforming minimum- weight path Py{a) first 
sweeps the subtree u W T{y \ v\). After that, for every node Uj, < i < fc, the 
path Py{a) sweeps the subtree Vi ttl T{C(vi) \ Vi+i) as follows: first, it sweeps a 
subtree Vi l±) T{Wi), finishing at Vi, and then, starting at Vi, it sweeps the subtree 
Vi l±l T{Wi), for some disjoint bipartitioning W, W Wi = C{vi) \ Vi+i. Finally, the 
path Py{a) sweeps the subtree T(a), finishing at a. 

The optimal choice of bipartitionings can be found as follows. We construct a 
weighted directed layered graph with a source vertex corresponding to node u = vo, 
a sink vertex corresponding to node v^. = a, and fc — 1 intermediate layers of 
vertices, each layer corresponding to a node Vi, < i < k. Each intermediate layer 
consists of at most 2*^"^ vertices, representing all different disjoint bipartitionings 
of the node set C{vi) \ Vi+i. The source and the sink vertices represent the trivial 
bipartitionings (l> \i) {V \ Vi) = V \ vi and C(a) W = C(a), respectively. Every 
consecutive pair of vertex layers (including the source and the sink vertices) are 
fully connected by forward arcs. In particular, the arc from a vertex representing 
the bipartitioning X ^ X in layer to the vertex representing the bipartitioning 
yi±)y in layer is given the weight ^{vi+\). It is easy to see that an optimal 
choice of bipartitioning corresponds to the minimum-weight path from the source 
to the sink in the layered graph. This minimum-weight path can be found by a 
standard dynamic programming algorithm (such as the Bellman-Ford algorithm, 
see e.g. [Gormen et al. 2001]) in time proportional to the number of arcs in the 
layered graph. 

Let W\ W Wi, . . . , Wk-i tbi Wk-i now denote the k — \ obtained optimal subtree 

bipartitionings. The k arcs of the corresponding source-to-sink shortest path deter- 
mine k edges (not necessarily consecutive) in the minimum-weight sweeping path 
Pyia). These edges are shown in Figure 3 by dotted lines. It now remains to apply 

the downsweep algorithm recursively in each of the subtrees u\±iT{V\vi), ?;il+lT(VKi), 
Vi^T{Wi), V2tiJT{W2), V2^T{W2), . . . , Vk-i^T{Wk-i), Vk-itiJT{Wk-i),T{a). 

Theorem 2.2. Given the output and the necessary intermediate values of the 
upsweep algorithm, the downsweep algorithm, computes the edges of the minimum- 
weight tree shortcutting in time and space (9(4'^n). 

Proof. The construction of the layered graph and the minimum-weight path 
computation runs in time 0(4'' /e), where k is the number of edges in the tree path 
w = t;o —> i"! —> V2 Wfc = a in the current level of recursion. Since the 

tree paths in different recursion levels arc edge-disjoint, the total number of edges 
in these paths is at most n. Therefore, the time complexity of the downsweep 
algorithm is 0{4'^n). 

By Theorem 2.1, the space complexity of the upsweep algorithm is 0{2'^n). In 
addition to the storage used internally by the upsweep algorithm, we also need to 
keep all the values DyT^{v). The number of possible pairs w, v is at most n (since 
node u is determined uniquely as the parent of v). For every such pair, we need to 
keep at most 4'' values, corresponding to 2'^ possible choices of each of V, W. The 
remaining space costs are negligible. Therefore, the total space complexity of the 
downsweep algorithm is 0{A'^n). □ 
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3. HEURISTIC IMPROVEMENTS 

Despite the guaranteed approximation ratio of the double-tree shortcutting and 
Christofides methods, neither has performed well in previous computational experi- 
ments (see [Johnson and McGeoch 1997; Reinelt 1994]). However, to our knowledge, 
none of these experiments explored the minimum-weight double-tree shortcutting 
approach. Instead, the double-tree shortcutting was performed in some subopti- 
mal, easily computable order, such as a depth-first tree traversal. We shall call this 
method depth-first double-tree shortcutting. 

In particular, [Reinelt 1994] compares 37 tour-constructing heuristics, including 
the depth-first doublc-trec algorithm and the Christofides algorithm, on a set of 
24 geometric instances from the TSPLIB database [Reinelt 1991]. Although most 
instances in this experiment are quite small (2000 or fewer points), they still allow 
us to make some qualitative judgement about the approximation quality of different 
heuristics. Depth- first double-tree shortcutting turns out to have the lowest quality 
of all 37 heuristics, while the quality of the Christofides algorithm is somewhat 
higher, but still far from the top. 

Intuitively, it is clear that the reason for the poor approximation quality of the 
two algorithms may be in the wrong choice of the shortcutting order, especially con- 
sidering that the overall number of alternative choices is typically exponential. This 
observation motivated us to implement the minimum- weight double-tree shortcut- 
ting algorithm from [Burkard et al. 1998] . It came as no surprise that this algorithm 
showed higher approximation quality than all the tour constructing heuristics in 
Reinelt's experiment. Unfortunately, Reinelt's experiment did not account for the 
running time of the algorithms under investigation. The theoretical time com- 
plexity of the previous minimum- weight doublc-trec algorithm from [Burkard ct al. 
1998] is 0{n^ -\- in practice, our implementation of this algorithm exhibited 

quadratic growth in running time on most instances. Both the theoretical and the 
practical running times were relatively high, which raised some justifiable doubts 
about the overall superiority of the method. 

As it was expected, the introduction of the new efficient minimum- weight double- 
tree algorithm described in Section 2 significantly improved the running time in our 
computational experiments. However, this improvement alone was not sufHcient 
for the algorithm to compete against the best existing tour-constructing heuristics. 
Therefore, wc introduced two additional heuristic improvements, one aimed at in- 
creasing the algorithm's speed, the other at improving its approximation quality. 

The first heuristic, aimed at speeding up the algorithm, is suggested by the well- 
known bounded neighbour lists [Johnson and McGeoch 2002, p. 408]. Given a tree, 
we define the tree distance between a pair of nodes a, b, as the number of edges 
on the unique path from o to 6 in the tree. Given a parameter k, the depth-k 
list of node u includes all nodes in the subtree T{u) with the tree distance from 
u not exceeding k. The suggested heuristic improvement is to limit the search 
across a subtree rooted at u in (l)-(2) to a depth-A; list of u for a suitably chosen 
value of k. Our experiments suggest that this approach improves the running time 
dramatically, without a significant negative effect on the approximation quality. 

The second heuristic, aimed at improving the algorithm's approximation quality, 
works by expanding the space of the tours searched, in the hope of finding a better 
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solution in the larger space. Let T be a (not necessarily minimum) spanning tree, 
and let A(T) be the set of all tours conforming to T, i.e. the exponential set of all 
tours considered by the doublc-trce algorithm. Our goal is to construct a new tree 
Ti, such that its node degrees are still bounded by a constant, but A(T) C A(ri). 
We refer to the new set of tours as an enlarged tour neighbourhood. 

Consider a node u in T, and suppose u has at least one child v which is not a 
leaf. We construct a new tree Ti from T by applying the degree-increasing operation, 
which makes node v a leaf, and redefines all children of v to be children of u. It is 
easy to check that any tour conforming to T also conforms to Ti. In particular, the 
nodes of T{v), which are consecutive in any conforming tour of T, are still allowed 
to be consecutive in any conforming tour of Ti. Therefore, A(T) C A(Ti). On the 
other hand, sequence w, u, v, where w is a child of v, is allowed by Ti but not by 
T. Therefore, A(T) C A(ri). 

Note that the degree-increasing operation cannot be performed partially: it would 
be wrong to reassign only some, instead of all, children of node u to a new parent. 
To illustrate this statement, suppose that v has two children wi and W2, which are 
both leaves. Let W2 be redefined as a new child of u. The sequence v,W2,wi is 
allowed by T but not by Ti, since it violates the requirement for v and W2 to be 
consecutive. Therefore, A(r) g A(Ti). 

We apply the degree-increasing heuristic as follows. Let D he a. global parame- 
ter, not necessarily related to the maximum node degree in the original tree. The 
degree-increasing operation is performed only if the resulting new degree of vertex 
u would not exceed D. Given a tree, the degree increasing operation is applied re- 
peatedly to construct a new tree, obtaining an enlarged tour neighbourhood. In our 
experiments, we used breadth-first application of the degree increasing operation 
as follows: 

Root the minimum spanning tree at a node of degree 1; 
Let r' denote the unique child of the root; 
Insert all children of r' into queue Q; 
while queue Q is not empty do 

extract node v from Q; 

insert all children of v into Q; 

if deg{parent{v)) -l-deg(t') < D then 

redefine all children of v to be children of parent{v) 

To incorporate the described heuristics, the minimum- weight double- tree algo- 
rithm from Section 2 was modified to take two parameters: the search depth k, and 
the degree limit D. We refer to the double- tree algorithm with fixed parameters k 
and D as a double-tree heuristic DTD,k- We use DT without subscripts to denote 
the original minimum-weight double-tree algorithm, equivalent to DTi^qq. 

4. COMPUTATIONAL EXPERIMENTS 

We compared experimentally the eSiciency of the original algorithm DT with the 
efficiency of double-tree heuristics DT^ k for two different search depths k = 16, 32, 
and for four different values for the degree limit D = 1 (no degree increasing 
operation applied), 3, 4, 5. The case D = 2 is essentially equivalent to D = 1, and 
therefore not considered. 
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5.97 


6.27 


6.43 


6.51 


6.47 




DT5,32 


5.62 


5.89 


5.93 


6.23 


6.38 


6.46 


6.43 





(a) Average excess over the Held-Karp bound (%) 



Size 


1000 


3162 


lOK 


31K 


lOOK 


316K 


IM 


3M 


DT 


0.18 


1.56 


15.85 


294..38 


3533 


51147 


156659 




DTi,i6 


0.04 


0.14 


0.47 


1.57 


5.60 


20.82 


101.09 


388.52 


DT3,i6 


0.10 


0.33 


1.12 


3.55 


11.90 


40.91 


138.41 


491.58 


DT3_32 


0.18 


0.69 


2.45 


7.56 


25.46 


82.99 


269.73 


935.55 


DT4,16 


0.23 


0.84 


2.78 


8.81 


29.02 


94.36 


307.31 




DT4.32 


0.45 


2.00 


6.93 


22.11 


74.70 


236.33 


744.50 




DT5,16 


0.62 


2.30 


7.79 


24.48 


81.35 


253.59 


807.74 




DT5,32 


1.11 


5.74 


20.73 


65.96 


224.34 


695.03 


2168.95 








(b) Average 


normalised 


running 


time (s) 







Table I: Results for DT and DT^i on uniform Euclidean distances 



The DIMACS Implementation Challenge [Johnson and McGeoch 2002] provided 
an excellent opportunity for testing and evaluating new approaches to the TSP. 
Website [DIMACS] , created to support the Challenge, contains a wide range of test 
instances and experimental data. In our computational experiments, we used uni- 
form random Euclidean instances with 1000 points (10 instances), 3162 points (five 
instances), 10000 points (three instances), 31623 and 100000 points (two instances 
of each size), 316228, 1000000, and 3168278 points (one instance of each size). 

For each heuristic, we consider both its approximation quality and running time. 
We say that one heuristic dominates another, if it is superior in both these respects. 
Following the approach of the DIMACS Challenge, approximation quality is mea- 
sured in terms of the approximate solution's excess over the Held-Karp bound 
(the solution to the standard linear programming relaxation of the TSP), and the 
running time in terms of the "normalised computation time" (sec [Johnson and Mc- 
Geoch 2002], [DIMACS] for details). The experimental results, presented in Table I, 
clearly indicate that nearly all considered heuristics^ (excluding DTi^ie) dominate 
plain DT. Moreover, all these heuristics (again excluding DTi^ig) dominate DT on 
each individual instance used in the experiment. 

For further comparison of the double-tree heuristics with existing tour-constructing 
heuristics, we chose DTi ie and DT5.16. 

The main part of our computational experiments consisted in comparing the 
double-tree heuristics against the most powerful existing tour-constructing heuris- 
tics. As a base for comparison, we chose the heuristics analysed in [Johnson and 



^Heuristic DTi_32 is omitted from Table I, since it does not give any noticeably better results 
compared to DTi^is. 
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Size 




1000 


3162 


lOK 


31K 


lOOK 


316K 


IM 


3M 


RA+ 




13.96 


15.25 


15.04 


15.49 


15.43 


15.42 


15.48 


15.47 


Chr-S 




11.48 


14.61 


14.81 


14.67 


14.70 


14.49 


14.59 


14.54 


r i 




12. .51 


12.17 


1:135 


1.3.11 


1.3. 39 


13. 13 


13. 17 




Sav 




11.38 


11.78 


11.82 


12.09 


12.14 


12.14 


12.14 


12.10 


ACh 




11.13 


11.00 


11.05 


11.39 


11.24 


11.19 


11.18 


11.11 


Chi-G 




9.80 


9.79 


9.81 


9.95 


9.85 


9.80 


9.79 


9.75 


Chr-HK 




7.55 


7.33 


7.30 


6.74 


6.86 


6.90 


6.79 




MTSl 




6.09 


8.09 


6.23 


6.33 


6.22 


6.20 






MTS3 




5.26 


5.80 


5.55 


5.69 


5.60 


5.60 






DTi,i6 




8.64 


9.24 


9.10 


9.43 


9.74 


9.66 


9.72 


9.66 


DT5,i6 




5.67 


5.91 


5.97 


6.27 


6.43 


6.51 


6.47 








(a) Average excess over 


the Held-Karp bound (%) 




Size 


1000 


3162 


lOK 


31K 


lOOK 


316K 


IM 


3M 


RA+ 


0.06 


0.23 


0.71 


1.9 


5.7 


13 


60 


222 


Chr-S 


0.06 


0.26 


1.00 


4.8 


21.3 


99 


469 


3636 


FI 


0.19 


0.76 


2.62 


9.3 


27.7 


65 


316 


1301 


Sav 


0.02 


0.08 


0.26 


0.8 


3.1 


21 


100 


386 


ACh 


0.03 


0.12 


0.44 


1.3 


3.8 


28 


134 


477 


Chr-G 


0.06 


0.27 


1.04 


5.1 


21.3 


121 


423 


3326 


Chr-HK 


1.00 


3.96 


14.73 


51.4 


247.2 


971 


3060 




MTSl 


0.37 


2.56 


17.21 


213.4 


1248 


11834 






MTS3 


0.46 


3.55 


24.65 


989.1 


2063 


21716 






DTi,i6 


0.04 


0.14 


0.47 


1.57 


5.60 


20.82 


101 


389 


DT5,i6 


0.62 


2., 30 


7.78 


24.48 


81.35 


254 


808 





(b) Average normaUsed running time (s) 



Table II: Comparison between established heuristics and DT-heuristics on uniform Euclidean 
instances 

Chr-S 

Sav 
VAch 

^ Chr-G 
"dTi.ib 

^DT3,i6 . Ch'-HK^ 

DT5,16 



2 4 6 8 10 12 14 16 
Average normalised running time (s) 

Fig. 4: Comparison between established heuristics and DT-heuristics on uniform Euclidean in- 
stances with 10000 points 
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McGeoch 2002], as well as two recent matching- based heuristics from [Kahng and 
Reda 2004]. The experiments were performed on a Sun Systems Enterprise Server 
E450, under SunOS 5.8, using the gcc 3.4.2 compiler. 

Table II shows the results of these experiments. Abbreviations in the table follow 
[Johnson and McGeoch 2002; Kahng and Reda 2004]: 

— RA+: Bentley's random augmented addition heuristic; 

— Chr-S: the Christofidcs heuristic with standard shortcut, implemented by John- 
son and McGeoch (JM); 

— FI: Bentley's farthest insertion heuristic; 

— Sav: saving heuristic, implemented by JM; 

— ACh: approximate Christofides heuristic, implemented by JM; 

— Chr-G: the Christofides heuristic with greedy shortcut, implemented by JM; 

— Chr-HK: the Christofides heuristic on Held-Karp trees instead of MST, imple- 
mented by Rohe; 

— MTSl, MTS3: "match twice and stitch" heuristics, implemented by Kahng and 
Reda. 

As seen from the table, the average approximation quality of DTi ig turns out to 
be higher than all classical heuristics considered in [Johnson and McGeoch 2002], 
except Chr-HK. Moreover, heuristic DTi ig dominates heuristics RA"'", Chr-S, 
FI, Chr-G. Heuristic DT5 ig dominates Chr-HK. Heuristic DT5 ig also compares 
very favourably with MTS heuristics, providing similar approximation quality at a 
small fraction of the running time. The above results show clearly that double-tree 
heuristics deserve a prominent place among the best tour-constructing heuristics 
for Euchdean TSP. 

The impressive success of double-tree heuristics must, however, be approached 
with some caution. Although the normalised time is an excellent tool for comparing 
results reported in different computational experiments, it is only an approximate 
estimate of the exact running time. According to [Johnson and McGeoch 2002, 
page 377], "[this] estimate is still typically within a factor of two of the correct 
time" . Therefore, as an alternative way of representing the results of computational 
experiments, we suggest a graph of the type shown in Figure 4, which compares 
the heuristics' average approximation quality and running time on random uniform 
instances with 10000 points. A normalised time t is represented by the interval 
[t/2,2t]. The relative position of heuristics in the comparison and the dominance 
relationships can be seen clearly from the graph. Results for other instance sizes 
and types are generally similar. 

Additional experimental results for clustered Euclidean instances are shown in 
Table III (with DTi^ig replaced by DT4^i6 to illustrate more clearly the overall 
advantage of DT- heuristics), and for TSPLIB instances in Table IV. 

While we have done our best to compare the existing and the proposed heuristics 
fairly, we recognise that our experiments are not, strictly speaking, a "blind test": 
we had the results of [Johnson and McGeoch 2002] in advance of implementing 
our method, and in particular of selecting the top DT-heuristics for comparison. 
However, we never consciously adapted our choices to the previous knowledge of 
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Size 


1000 


3162 


lOK 


31K 


lOOK 


316K 


RA+ 


12.84 


13.88 


16.08 


15.59 


16.22 


16.33 


Chr-S 


12.03 


12.79 


13.08 


13.47 


13.50 


13.45 


I'l 


9.90 


11.85 


12.82 


1. ■-!.:-! 7 


i:i.9(i 


13.92 


Sav 


13.51 


15.97 


17.21 


17.93 


18.20 


18.50 


ACh 


10.21 


11.01 


11.47 


11.78 


12.00 


11.81 


Chr-G 


8.08 


9.01 


9.21 


9.47 


9.55 


9.55 


Chr-HK 


7.27 


7.78 


8.37 


8.42 


8.46 


8.56 


MTSl 


8.90 


9.96 


11.97 


11.61 


9.45 




MTS3 


8.52 


9.5 


10.11 


9.72 


9.46 




DT4,16 


6.37 


8.24 


8.79 


9.40 


9.38 


9.39 


DTs.ie 


5.72 


7.17 


7.92 


8.32 


8.46 


8.42 


(a) Average excess over the Held-Karp bound (%) 


Size 


1000 


3162 


lOK 


31K 


lOOK 


316K 


RA+ 


0.1 


0.2 


0.7 


1.9 


5.5 


12.7 


Chr-S 


0.2 


0.8 


3.2 


11.0 


37.8 


152.8 


FI 


0.2 


0.8 


2.9 


9.9 


30.2 


70.6 


Sav 


0.0 


0.1 


0.3 


0.9 


3.4 


22.8 


ACh 


0.0 


0.2 


0.8 


2.1 


6.4 


54.2 


Chr-G 


0.2 


0.8 


3.2 


11.0 


37.8 


152.2 


Chr-HK 


0.9 


3.3 


11.6 


40.9 


197.0 


715.1 


MTSl 


0.78 


4.19 


45.09 


276 


1798 




MTS3 


0.84 


4.76 


49.04 


337 


2213 




DT4,16 


0.2 


0.87 


3.16 


9.55 


34.43 


120.3 


DT5,16 


1.12 


4.85 


16.08 


53.35 


174 


569 



(b) Average normalised running time (s) 



Table III; Comparison between established heuristics and DT-heuristics on clustered Euclidean 
instances 

[Johnson and McGeoch 2002], and we believe that any subconscious effect of this 
previous knowledge on our experimental setup is negligible. 

5. CONCLUSIONS AND OPEN PROBLEMS 

In this paper, we have presented an improved algorithm for finding the minimum- 
weight double-tree shortcutting approximation for Metric TSP. We challenged our- 
selves to make the algorithm as efficient as possible. The improvement in time 
complexity from 0{n^ + to 0(4''n^) (which implies O(n^) for the Euclidean 

TSP) placed the minimum- weight double-tree shortcutting method as a peer in the 
set of the most powerful tour-constructing heuristics. It is known that most such 
heuristics have theoretical time complexity O(n^), and in practice often exhibit 
near-linear running time. The minimum-weight double-tree method now also fits 
this pattern. 

While we have not been using the language of parameterised complexity [Downey 
and Fellows 1998], wc (and the previous work [Burkard et al. 1998]) have in fact 
demonstrated that the problem of finding the minimum- weight double-tree tour for 
Metric TSP is fixed-parameter tractable (where the maximum degree of the MST 
is the relevant parameter). It would be interesting to see if this connection with 
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Size 


1000 


3162 


lOK 


31K 


lOOK 


RA+ 


17.46 


16.28 


17.78 


19.88 


17.39 


Chr-S 


13.36 


11.17 


13.41 


16. .50 


15.16 


¥1 


15. T)!) 


11.28 


13.20 


17.78 


15.32 


Sav 


11.96 


12.14 


10.85 


10.87 


19.96 


ACh 


9.64 


10.50 


10.22 


11.83 


11.52 


Chr-G 


8.72 


9.41 


8.86 


9.62 


9.50 


Chr-HK 


7.38 


7.12 


7.50 


6.90 


7.42 


MTSl 


7.0 


6.9 


5.1 


4.7 


4.1 


MTS3 


6.2 


5.1 


4.0 


2.9 


2.7 


DTi,i6 


6.36 


5.99 


8.09 


9.99 


10.02 


DTs.ie 


6.13 


5.58 


7.65 


8.98 


9.30 


(a) Average excess 


over the Held-Karp bound (%) 


Size 


1000 


3162 


lOK 


31K 


lOOK 


RA+ 


0.1 


0.2 


0.8 


2.2 


5.6 


Chr-S 


0.1 


0.2 


1.8 


3.9 


31.8 


FI 


0.2 


0.8 


3.1 


9.8 


26.4 


Sav 


0.0 


0.1 


0.3 


0.6 


1.4 


ACh 


0.0 


0.1 


0.5 


1.5 


3.9 


Chr-G 


0.1 


0.2 


1.8 


3.8 


29.5 


Chr-HK 


0.7 


2.2 


9.7 


50.1 


177.9 


MTSl 




1.5 


34.4 


107.3 


620.0 


MTS3 




2.1 


42.4 


135.4 


1045.3 


DTi,i6 


0.3 


0.9 


4.1 


18.4 


49.3 


DTs.ie 


0.6 


2.1 


11.0 


57.1 


115.1 



(b) Average normalised running time (s) 



Table IV: Comparison between established heuristics and DT-heuristics on geometric instances 
from TSPLIB: prl002, pcbll73, rll304, nrwl379 (size 1000), pr2392, pcb3038, fnll4461 (size 
3162), pla7397, brdl4051 (size lOK), pla33810 (size 31K), pla859000 (size lOOK). 

parameter ised complexity theory can be extended further, e.g. by using any of the 
estabUshed tccliniques for designing fixed-parameter tractable algorithms. 

Our results should be regarded only as a first step in exploring new opportunities. 
Particularly, the minimum spanning tree is not the only possible choice of the initial 
tree. Instead, one can choose from a variety of trees, e.g. Held and Karp (l-)trees, 
approximations to Steiner trees, spanning trees of Delaunay graphs, etc. This 
variety of choices merits a further detailed exploration. 

It is "well-known that when the initial tree is a path, the resulting double-tree 
tour neighborhood is the set of all pyramidal tours [Burkard et al. 1998]. In this 
case, a dozen of conditions on the distance matrix are known (see e.g. [Burkard 
et al. 1998]), which guarantee that the tour neighbourhood contains the absolute 
minimum- weight tour. It may be possible to generalise this approach by identi- 
fying new special types of trees and conditions on the distance matrices, which 
would guarantee that the minimum- weight double-tree algorithm finds an absolute 
minimum-weight tour. For more results on polynomial solvability of TSP with spe- 
cial conditions imposed on the distance matrix, see [Burkard et al. 1998; Deineko 
et al. 2006]. 

ACM Journal Name, Vol. V, No. N, Month 20YY. 



Fast minimum-weight double-tree shortcutting for IVIetric TSP 



15 



The minimum-weight shortcutting problem for the Christofides graph remains 
NP-hard even in the planar Euclidean metric. However, our algorithm turns out 
to be applicable also to this problem on certain classes of instances. It can be 
shown that if the Christofides graph is a cactus (i.e. all its cycles are pairwise edge- 
disjoint), then the set of all its shortcuttings is a subset of the set of all double- tree 
shortcuttings. Therefore, our algorithm, as well as the algorithm of [Burkard ct al. 
1998], can be used to find efficiently the minimum- weight shortcutting when the 
Christofides graph is a cactus. In particular, such a shortcutting can be found in 
polynomial time in the planar Euclidean metric. 

Our efforts invested into theoretical improvements of the algorithm, supported by 
a couple of additional heuristic improvements, have borne the fruit: computational 
experiments with the minimum- weight double-tree algorithm show that it becomes 
one of the best known tour constructing heuristics. It appears that the double-tree 
method is also well suited for local search improvements based of transformations of 
trees and searching the corresponding tour neighborhoods. One can easily imagine 
various tree transformation techniques that could make our method even more 
powerful. 
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